3D simulations of Einstein's equations: symmetric hyperbolicity, live gauges and 

dynamic control of the constraints 
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We present three-dimensional simulations of Einstein equations implementing a symmetric hy- 
perbolic system of equations with dynamical lapse. The numerical implementation makes use of 
techniques that guarantee linear numerical stability for the associated initial-boundary value prob- 
lem. The code is first tested with a gauge wave solution, where rather larger amplitudes and 
for significantly longer times are obtained with respect to other state of the art implementations. 
Additionally, by minimizing a suitably defined energy for the constraints in terms of free constraint- 
functions in the formulation one can dynamically single out preferred values of these functions for 
the problem at hand. We apply the technique to fully three-dimensional simulations of a station- 
ary black hole spacetime with excision of the singularity, considerably extending the lifetime of the 
simulations. 



I. INTRODUCTION 

The construction of an accurate and stable numerical 
implementation of Einstein equations in settings such as 
three dimensional binary black hole collisions represents 
a major challenge. Indeed, although significant progress 
has been achieved in the last few years (see jy and ref- 
erences therein), the goal remains elusive. In practice, 
numerical or continuum instabilities often arise, render- 
ing particular simulations of little use after some time, 
or at high enough resolutions. Due to the complicated 
set of equations one is dealing with, coupled to the often 
scarce computational resources, tracking down the source 
of problems is usually difficult. Faced with this situation, 
a possible strategy is to make use of techniques that sys- 
tematically control different aspects of the problem under 
consideration. One possibility for such a strategy can be 
achieved by proceeding along the following lines: 

First, as has been emphasized in a number of works 
, by choosing a well posed initial-boundary value prob- 
lem (IB VP), which is a necessary condition for a numer- 
ically stable implementation (see, for example, 0,3). A 
symmetric hyperbolic system with maximally dissipative 
boundary conditions is an example of a system that de- 
fines such a well posed problem. Furthermore, the bound- 
ary conditions should not only define a well posed initial 
value problem but also must conform to the physical sit- 
uation in mind. For instance, they should preserve the 
constraints and have minimal spurious influence on the 
solution. 

Second, by translating the previous analytically- 
framed considerations into the numerical arena. That 
is, by constructing a numerically stable scheme for the 
IBVP under consideration. One way of doing so is by 
constructing difference operators and imposing the dis- 
crete boundary conditions in a way such that the steps 
followed at the continuum to show well_posedness can be 
reproduced at the discrete level 



Third, by implementing the equations so that spurious 
growth in time of the solution is removed or minimized. 
A numerically stable implementation need not free the 
simulation from errors that at fixed resolution grow fast 
in time (although at fixed time they would go away with 
resolution). There are a number of possible alternatives 
for this effect to be minimized. 

One option is to achieve semi-discrete or discrete strict 
stability, i.e. given a sharp energy estimate at the contin- 
uum, to discretize in a way such that the estimate also 
holds at the semi-discrete or discrete levels 0, ■ ^his 
way, growth in the numerical solution that is not called 
for by the continuum system is ruled out. In this strat- 
egy one must count with a sharp energy estimate for the 
problem at hand, such as those of [8ll3.ITol|. 

An alternative and/or complementary way is to con- 
sider the addition of a small amount of artificial dissipa- 
tion (in a way that does not spoil the available discrete 
energy estimates and numerical stability). Sometimes 
this is enough to partially or completely rule out errors 
growing fast in time, and in addition it helps to control 
high frequency modes. 

There are other cases, in which a sharp estimate is 
not available, and for which the addition of some artifi- 
cial dissipation does not rule out undesired growth in the 
solution either. This happens quite often in evolutions 
of Einstein's equations in the strong field regime. One 
possibility in such a case is to realize that even though 
a sharp energy estimate for the main evolution system 
might not be available, an ideal one for the subsidiary 
system that describes how the constraints propagate is 
trivially obtained Namely, at the continuum one 

would like the constraints when perturbed to (for exam- 
ple) remain constant as a function of time, or decay to 
zero [T^ . Similarly, at the discrete level one may want 
constraint violations to remain close to their initial, dis- 
cretization value. Once some desired estimate for the 
constraints is chosen, one can enforce it in a number of 
ways. One of them is to dynamically redefine the equa- 
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tions during evolution off of the constraints surface ■ 
In this paper we present results obtained with a fully 
non-linear code that evolves Einstein's equations in a 
three-dimensional setting, analyzing implementations of 
techniques that ensure some of the desired properties just 
discussed. The particular formulation of the equations 
that we use is a symmetric hyperbolic one with a dy- 
namical gauge condition (more precisely, a slight gener- 
alization of the Bona-Masso slicing conditions) presented 
in Ref. [l^ , and summarized in Section ^] In Section 
mil details on how free constraint-functions in the for- 
mulation can be dynamically adjusted to ensure mini- 
mal growth of some energy, or norm, associated with 
constraints, is discussed in the context of the symmet- 
ric hyperbolic formulation here used. Section llVI brieflv 
summarizes the numerical techniques used in this pa- 
per, already presented in Ref. [g, and the details of 
the test-beds here studied. One of these test-beds is the 
study of a periodic gauge wave, presented in Section Ivl 
There we show that the use of a symmetric hyperbolic 
formulation and a small amount of artificial dissipation 
suffices to evolve this solution with rather large ampli- 
tudes and for long times. The other test-bed that we 
study is a non-spinning, stationary black hole with exci- 
sion of the singularity and dynamic minimization of the 
constraints' growth. A detailed analysis is presented in 
Section IVll discussing several issues relevant to the con- 
straint minimization technique and the results of fully 
three-dimensional simulations whose lifetime is consider- 
ably extended by making use of this technique. Section 
IVIII summarizes and discusses the main lessons of this 
work and points out possible extensions of it. 



II. THE SYMMETRIC HYPERBOLIC 
FORMULATION USED 

In this paper we use the symmetric hyperbolic formu- 
lation of the Einstein equations admitting a dynamical 
lapse introduced in p^ . This system has thirty- four 
variables, including: the three metric, gij, the extrin- 
sic curvature, Kij, and the lapse, N. Further, variables 
dkij and Ai are constructed from the spatial derivatives 
of gij and N, respectively, and introduced as indepen- 
dent variables to make the system first order in space. 
When all constraints are satisfied these variables satisfy 
dkij = dkSij and Ai = N~^diN. The evolution equations 
in this formulation are: 



dogij 
doKij 

dodkij 
doN 



= -2dkK,j - 2AkKij + 77(0;^) gk{iCjf 
= -F{N,K,x^') + Six") 



(1) 



(2) 

(3) 
(4) 



ON 

1 dF{N,K,x'') 



N 



dx^ 



N dK 

+ ^{xnC^, 



(5) 



where we define do = N ^ {dt — £ p) ■ The Ricci tensor 
in Eq. (|2Jl is written as 



Ri 



9°"^ {-dadbij + dad(ij)i, + d(^.id\ab\]) - d{idj)ah) 



d,^^d 



'jab 



{dk - 2bk)r% - T)^T\, , 



where bj = dkijg , dk = dki^g^^ , and 



T% = ig'^' (2d 



Finally, the second order derivatives of N that appear in 
Eq. 101 are calculated as 



TV 



r ijAk + AiAj 



The slicing condition, eq. Q), contains two functions, 
F{N,K,x^) and S{x^). The function F may be any 
function of the lapse, the trace of the extrinsic curva- 
ture, K = g^^Kij, and the spacetime coordinates, with 
the condition that dF/dK > 0. The function 5 is a 
gauge source function, and is specified a priori but in 
an arbitrary way as a function only of the spacetime co- 
ordinates. This slicing is a generalization of the Bona- 
Masso slicing conditions, obtained by setting 5 = and 
F{N,K,x^) = f{N)K. Moreover, choosing 5 = and 
f = N gives the harmonic time slicing condition, or a 
generalized harmonic condition if 7^ 0. The time har- 
monic slicing condition, or its generalized form, is the 
choice used for all runs in this paper. Finally, the shift, 
(3^{x^), must also be specified a priori as an arbitrary 
function of spacetime. 

The Einstein equations are a constrained system, and 
the evolution equations here considered are not only sub- 
ject to the physical constraints, the Hamiltonian and 
momentum ones, but also non-physical constraints that 
come from introducing first-order variables. The Hamil- 
tonian constraint is 



C=iR-KabK'''' + K^)/2, 



(6) 



where R = g^-' Rij and the Ricci tensor given by Eq 
The momentum constraints, Ci = V^Kai — ViK, are 



a = 9"' (daKM - mab) + 2 (rf' - 26')iffe. + ^'Kab ■ 

(7) 

Finally, the non-physical constraints, Ca^ , Ckij and Cikij , 
are defined as 



Ckij 
Clkij 



A, - N-^m , 

d[idk]ij. 
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The Einstein equations resulting from the 3+1 ADM 
decomposition are only weakly hyperbolic. However, it 
is possible to manipulate the principal part of the equa- 
tions by adding the constraints in specific combinations 
to the evolution equations in order to obtain a strongly 
or symmetric hyperbolic system of equations |14| . In 
the system here considered, the constraints are added 
to the RHS of Eqs. (^3)), and the spacetime constraint- 
functions {jXtViX^^} Eire introduced as multiplicative 
factors to the constraints. Requiring the evolution sys- 
tem to be symmetric hyperbolic imposes algebraic con- 
ditions on these factors, as discussed below, and they are 
not treated as completely independent. Typically these 
factors are taken to be constant parameters, however this 
restriction is actually not needed for strong or symmetric 
hyperbolicity of the system to hold. Moreover, we wish 
to exploit some freedom in choosing these constraint- 
functions to minimize the effect of constraint violating 
modes that may appear in the solution llj . Thus, we 
choose the factors to be functions of time but constant 
in space (future work will concentrate in allowing for 
space dependence): {7(0: C(*)7 ^(0: There- 
fore, here we will refer to these factors as constraint- 
functions rather than parameters. 

The characteristic speeds of the system are (3^ni,±N+ 
P'n^, ±Ny/X; + fihi,, with 

Ai = 2cre//, 

A2 = l + x-^(l + C)^ + 7(2-r? + 2x), 

A3 = '\x-\{K^i)v-\i- 

and (Jeff = {dF/dK)/{2N). 

There arc two strongly hyperbolic multiple constraint- 
function families with A2 = 1 and A3 = 1, which cor- 
respond to propagation speeds along the light cone and 
the hypersurface normal. One such family has three free 
constraint-functions {7(t), C,{t),ri{ty\: 



1 ^ 



X 



(1 + 0^-27(2-7?) 
2(1 + 27) 



C = -^X-^(3C-l)?/-2, 



A second family has only two free constraint-functions 
{C(i),xW}: 



1 . ^.111 



In this paper we set C = —1 to simplify the calcu- 
lation of the characteristic variables needed to impose 
maximally dissipative boundary conditions. This gives 
a symmetric hyperbolic system with one free constraint- 



function x{i) 



Single constraint-function system < 



r 7 = 


1 

2 


c = 


-1 




2 


e = 


X 

2 


I X ^ 






(8) 



and another symmetric system with two varying 
constraint-functions {r]{t),"f{t) ^ —1/2}: 



Two constraint-function system < 



r c = 

X = 

i = 

7 ^ 



-1 

1+27 



(9) 

One can show that these families not only are strongly 
hyperbolic as shown in Hal but, as mentioned above, 
also symmetric hyperbolic 



III. DYNAMIC CONTROL OF THE 
CONSTRAINTS' GROWTH. 

The formulation of the Einstein equations summarized 
in the previous section is made symmetric hyperbolic by 
adding constraints to the evolution equations multiplied 
by the time-dependent constraint-functions. Requiring 
that the propagation speeds be along the light cones or 
t =constant hypersurfaces normal, and a further simpli- 
fication obtained by setting C = — 1, results in two fam- 
ilies of equations. The first has a single free constraint- 
function, X, and the second has two constraint-functions, 
{7, 77}. Several papers have been presented showing that 
the long-term stability of 3D numerical simulations is ex- 
tremely sensitive to these constraint-functions, even if 
symmetric hyperbolicity is guaranteed independently of 
the values these constraint-functions take |lfl| . Recently, 
in ^1 a method to dynamically choose these constraint- 
functions by minimizing the constraint growth during the 
evolution has been presented. We here include a brief 
summary of this method, and discuss its particular ap- 
plication to the black hole runs presented later in this 
paper. 

Consider a system of hyperbolic equations with con- 
straint terms, Cc, written schematically 

Ua = '^A''{u,t,x)dbUa+Ba{u,t,x)V^pLacCc{u,djU) , 
b c 

(10) 

where Uq, Ba and Cc are vector valued functions, and /Iqc 
is a matrix (generally not square) that is a function of 
the spacetime. (Note that in this section Cc represents 
a vector function of general constraint variables, and not 
specifically the momentum constraint introduced in the 
previous section.) Here the indices {a,6, c} range over 
the size of each vector or matrix function, while the in- 
dices {i,i,k} will label points on a discrete grid. We 
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then define an energy or norm of the discrete constraint 
variables, e.g., 



N{t) 



1 



(11) 



where nx,riy,nz are the number of points in the x,y,z 
directions, and where we have omitted the grid indices 
k} to simplify the notation. The time derivative of 
the norm can be calculated using Eq. (|10|l 



rhojn 



Tr(/il^) 



(12) 



and therefore can be known in closed form provided the 
matrix valued sums 



Ca 



EE' 

C 

EE 

ijk a 

dCg 

duh 



dCg ^ 
Bb 



dCg 

dDkUa 



(14) 



are computed during evolution; where Di is the discrete 
derivative approximation to di. We then use the de- 
pendence of the energy growth on the free constraint- 
functions to achieve some desired behavior for the con- 
straints, i.e., solving Eq. (|12|l for jiac- For example, if we 
choose [33 

jif = -aAf, a > 0, (15) 

any violation of the constraints will decay exponentially 

M{t + M)=N{t)e-''^' . (16) 

As discussed in i one good option among many others 
seems to be choosing a tolerance T value for the norm of 
the constraints that is close to its initial, discretization 
value, and solving for ^ac such that the constraints decay 
to this tolerance value after a given relaxation time. This 
can be done by adopting an a such that after time riaAt 
the constraints have the value T. Replacing M{t + A<) 
by T in equation l|lt)|) and solving for a gives 



a{t) 



1 



If one then solves 

Af = -aAf = + trace(^ x X^) 



(17) 



(18) 



for /X, with a given by Eg. ljlTII . the value of the norm 
Af{t + na^t) should be T, independent of its initial value. 



In principle, Ea. (|18|l has non-unique solutions, since 
the equation is scalar and /ioc is a matrix. As discussed 
in Section fVI Bl this non-uniqueness is sometimes crucial 
for making this a practical method for controlling growth 
in the constraints. 

We also note that the technique discussed in this sec- 
tion can be implemented without affecting symmetric hy- 
perbolicity llj. 

Finally, we now describe how Af is calculated for the 
symmetric hyperbolic families of the Einstein equations 
used in this paper. The time derivative of the energy for 
the constraints is 

Af = + + + li^ + ii^ + rii'^ (19) 

For the single- function family, Eq. ||SJ), 



(13) That is. 



(20) 



with 



•j-hom j-hom ^^')' J^^^ -f\- 

2 

The evaluation of TV" as a function of x is a two-step 
process. In order to compute the quantities 1'*°™ and 

, so as to obtain the dependence of the time derivative 
of the energy in terms of x, Eq. H2U|) , two evaluations of N 
are required to extract the individual contributions. For 
example, the homogeneous term is obtained performing 
a set of evaluations with x = 0, 



rhom 



: ^(X = 0) 



Once this term is known, is obtained doing a set of 
evaluations with an arbitrary but non- vanishing x — XO: 



rhom 



Xo 



This involves at each step an evaluation of the right hand 
side of the evolution equations and of the spatial deriva- 
tives of such right hand side, and evaluation of the deriva- 
tives of the constraints with respect to the main variables 
and with respect to their spatial derivatives (all of this 
at each grid point) The constraints, and therefore their 
derivatives, do not depend on the constraint-functions. 
Therefore they need to be computed only once at any 
given time. 

Similarly for the two- function family l^: the time 
derivative of the norm of the constraints is 



Af- 



rhom 



Xl^ + + ryr' 



(21) 
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with [c.f. Eq. ®] 



Numerical method 



and 



fhom 



7(2 - V) 
1 + 27 



— Tx _ It" 
2 ' 



At any given time, four sets of evaluations are needed 
to numerically compute the quantities , X'^ , X''' , X'' 
in Eq. fill) . For example, as in the single constraint- 
function case, the homogeneous term is obtained through 
a set of evaluations with 77 = = 7, 

X'""" = iV,(77 = 0,7 = 0). 

Once this is known, X'' is obtained doing a set of evalu- 
ations with an arbitrary but non- vanishing rj = rjo, and 
7 = 0, 



X" = 



(22) 



Two more sets of evaluations are needed in order to con- 
struct X'' and X"*": given 71 and 72 arbitrary but different, 
71 7^ 72 , it is straightforward to see that 

j-y ^ (72 + 27i72)jVc(7i) - (72 + 27i72)iVc(7i) ^^^^ 



1^ 



27172(71 - 72) 
(l-H27i)(l-|- 272) 



47172(71 - 72) 



72A^c(7i) -7i^c(72) (24) 



The quantities X'', X^, thus obtained are indepen- 
dent of what values f?o, 7i, 72 are used in Eqs. (|22l23lEl|l . 
We make use of this fact to perform a non-trivial test of 
self-consistency in our simulations. Namely, during evo- 
lution we construct these quantities X**, X^, X^ using, 
at each time step, several different values of 770,71,72, 
and we check that the result is, indeed, independent of 
that choice. We proceed in a similar way with the single- 
constraint-function family. 



IV. NUMERICAL METHODS AND TEST 
PROBLEMS 



We use numerical techniques based on the energy 
method for hyperbolic equations . This method allows 
one to identify numerically stable discretizations by con- 
struction for initial-boundary value problems for linear, 
symmetric hyperbolic systems. While we focus here on 
nonlinear Einstein equations, we note that some numeri- 
cal instabilities in the Einstein system are also observed 
in the linear regime. Methods that are known to be nu- 
merically stable for the linearized Einstein equations thus 
function both as a foundation and guide for moving to the 
nonlinear problem. Our numerical scheme uses second 
order spatial difference operators that satisfy summation 
by parts, as well as an extension of the standard Kreiss- 
Oliger dissipation operator that takes into account the 
presence of (inner and outer) boundaries. This operator, 
which we call Qd^ is added to the right-hand side of the 
evolution equations, w = (. . .) -t-Qd, with a free (non neg- 
ative) multiplicative parameter a. Paper describes in 
detail this operator; here to fix ideas we include it for the 
non-excision case. It is, on each direction, 

Qd'^a = — 2(TAa;D^uo, 

Qdu, = -a{Axf{D+D^fu,, for ^ = 2, . . . , iV - 2 

QdUN-i = -<tAx{D^_ -2D+D_)uN-i, 

QdUN = —2aAxD^UN. (25) 

Maximally dissipative boundary conditions are im- 
posed numerically through projections that are orthog- 
onal in the linear case. We use third order Rungc- 
Kutta to integrate the equations in time. The compu- 
tational domain consists of an uniform Cartesian grid 
(Ax = Ay = Az = A). The black hole simulations 
employ a cubical inner boundary to excise the singular- 
ity from the computational domain. Unless otherwise 
stated, the simulations presented throughout the paper 
use a dissipative parameter a — 0.03 and Courant factor 
A — 0.5. This choice for a is motivated by the fact that, 
at least for the second order wave equation written in 
first order form, it gives the maximum Courant factor al- 
lowed by von-Neumann stability (see 0]). For additional 
information on our numerical scheme see 

Boundary conditions are specified via maximally dis- 
sipative boundary conditions where needed. These are 
introduced by finding all the incoming modes and setting 
the time derivative of these to zero. Maximally dissipa- 
tive boundary conditions can be written in the form 



In this section we introduce the numerical methods 
that we use, and then outline two physical problems, the 
gauge wave and a Schwarzschild black hole, that we will 
analyze in this paper. These spacetimes will be used 
in our numerical implementation of Einstein equations 
written in the first order form detailed in section II - 
equations through 



Lu +g{t,x^). 



(26) 



where L is "sufficiently small," such that an energy es- 
timate for the IBVP can be derived for symmetric hy- 
perbolic systems in more than one dimension (see, e.g., 
|17|.'). The function (7 is an a priori but arbitrary func- 
tion of the spacetime coordinates in the boundary and 
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time. The combination of maximally dissipative bound- 
ary conditions and a symmetric hyperbolic evolution sys- 
tem define a well-posed initial-boundary value problem. 
However, these boundary conditions in general are not 
constraint-preserving. For the evolutions presented here, 
setting the time derivative of all the incoming modes to 
zero is consistent with preserving the constraints at the 
boundary because the exact solution is known. However, 
this will not hold in general and future work will concen- 
trate on evolutions with constraint-preserving boundary 
conditions. 



B. Gauge waves 

We first test our numerical method by studying a gauge 
wave defined by 



(27) 



which corresponds to a coordinate transformation in the 
(x, t) plane, of flat spacetime. The analytic solution for 
the gauge wave is obtained by setting the gauge source 
function to zero, S{x^') = 0, in Eq. ®: -N^K = dtN 
and the shift — 0. We adopt periodic boundary con- 
ditions to simplify the analysis by eliminating possible 
boundary effects. 



C. Black Holes 

We then examine in some detail tests with 
Schwarzschild spacetime, which describes a static non- 
spinning black hole. The singularity inside the black hole 
is excised, which restricts the possible slicings (surfaces 
of constant time) that we consider to those that smoothly 
penetrate the black hole horizon, such that we can place 
an inner boundary inside of the horizon. In the present 
work we consider a Schwarzschild black hole in Kerr- 
Schild (KS) — or ingoing Eddington-Finkelstein — coor- 
dinates. In these coordinates, the metric of the spacetime 
is given by the line clement: 

ds^ = -N^df + grj{dx' + p'dt){dx^ + (3Ht) , (28) 



where. 



N 



r + 2m 



1/2 



P = — X 



r + 2m 

2m x^x^ 



(29) 
(30) 
(31) 



The gauge source function S is read off from the exact 
solution. 



V. GAUGE WAVE SIMULATIONS 

The gauge wave is a simple, non-trivial numerical test 
problem, as it is free of boundaries, the amplitude of 
the fields can be controlled by a single parameter, and 
it does not lead to any singularity. This solution is used 
to compare the performance of different implementations 
of Einstein's equations . Despite its simphcity at the 
analytical level, this test illustrates the challenges associ- 
ated with the numerical implementation of Einstein equa- 
tions 2^ In particular, it is often observed 
that for amplitudes A > 0.01, the numerical solutions 
display exponential growth and loss of convergence. 

Since in the cases discussed in this paper the analytical 
solution is known, the convergence factor can be defined 
as 



C{F) = \\F^-F, 



analytical \ \2/ 



I-Fa/2 - Fa. 



analytical \ \2 : 

(32) 

with F the variable under consideration, F/x and 
Fanaiyticai its numerical (at resolution A) and analyti- 
cal solutions, respectively. 

In the present context, the gauge wave is used for two 
purposes. First, it allows us to probe the stability of our 
numerical method, showing its advantages and present 
limitations. Second, it sheds light on possible sources 
of instabilities or spurious growth often encountered in 
these tests. In particular, we are able to establish that 
constraint violations, if any, grow quite slowly, allow- 
ing one to accurately follow the system for thousands 
of crossing times, even for large values of the amplitude 
parameter A. We have ran our simulations for a wide 
range of values for A, observing qualitatively the same 
results. 

Note that one can readily show that the constraints, 
C, Ci, and Cijki, for the gauge wave in this formula- 
tion are satisfied to round-off level if the variables depend 
solely on {t,x). That C, Ci are satisfied follows directly 
from the fact that — {f{x) — l)(5^<5j -I- r]ij, Kij = 
g{x)6lSi, d,jk = h{x)6l6iS^; with {g{x), f{x),h{x)} &t- 
bitrary functions. Note that: (i) at any t = const hyper- 
surface , gij describes a flat metric and so Rij = analyti- 
cally. Furthermore, at the numerical level all terms in 
cancel, including those describing the truncation errors in 
the derivatives, (ii) — KijK^^ and Kij —gijK are zero 
algebraically and so, coupled to (i), C, Ci are satisfied nu- 
merically to round-off level, (iii) Cijki is defined by the 
conmutator of two derivatives; since the only non-trivial 
components are in the x direction, the only surviving one 
is Cxxxx = d[xdx]xx which is trivially satisfied. Conse- 
quently, the value of the constraint-function multiplying 
them in the right hand side of the evolution equations do 
not play a significant role. This observation holds for the 
variables depending solely on {t, x), as mentioned, this is 
indeed the case throughout our runs. Therefore, for these 
tests we adopt the two-constraint-function family for- 
mulation with fixed values for the constraint-functions: 
7y = 7 = 0. 
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We concentrate on two non-linear cases with relatively 
large amplitudes: A = 0.5 and A = 1. For example, when 
A = 0.5, gxx ranges over the interval [0.6, 1.7], and over 
[0.37,2.72] when A — 1. In 7], an independent analysis 
and code for the linear equations around the gauge wave 
was presented, illustrating the tendency of the numerical 
solution to grow exponentially unless some small amount 
of dissipation is added to the RHS. 

The computational domain is chosen to be [—1,1], 
and is represented on a uniform grid with spacing A = 
2/{N — 1). The gauge wave is really a one dimensional 
problem, with a non-trivial dependence only on {x,t). 
The solution is constant in the (y, z) coordinates, thus 
only a few points are need to represent the field in these 
directions. We therefore use = {80p + 1) points with 
p = 1,2,4 in the x direction, and five points in each of 
the (y, z) directions (though tests were performed with 
uniform grids, equally spaced in all directions obtaining 
exactly the same results as expected). In the gauge wave 
tests, we chose a Courant factor of 1/4 so as to make a 
more direct comparison with similar tests presented in 
the literature (see for instance |1^ ilE 113 ) • 



1. Amplitude A — 0.5 

Figures El and 13 illustrate our results for the gauge 
wave with amplitude A — 0.5. Figure ^ shows for three 
different resolutions the logarithm of the energy for the 
constraints as a function of time, for 1250 crossing times. 
The lowest resolution shows a slowly growing behavior, 
but this is considerably diminished as better refined ones 
are considered. Note that even for the coarsest resolu- 
tion, the energy remains smaller than 10~^. 

In this type of simulation, phase differences between 
the numerical and exact solutions typically cause the er- 
ror to vary in an oscillatory way: the error going back 
to a small value after some time, when the numerical so- 
lution achieves a phase difference of 2tt relative to the 
exact one Q. To give a phase-independent indication of 
the errors with respect to the exact solution, in Figure 
121 we show the relative error in the maximum attained 
by gxx, compared to the exact value. A slow growth is 
observed. In particular, the error in the maximum of the 
function is ~ 0.07 even after 1250 crossing times for the 
finest discretization. Fig. (21 also explicitly shows that the 
amount of dissipation is very small, as the amplitude is 
not damped even in very long runs. The crossing of lines 
in Fig. 121 rfoes not imply lack of convergence, since one is 
not necessarily comparing the field at the same points. 
(For an explicit convergence illustration of the solution it- 
self, we present some results for the more stringent A = I 
case). Figure n explicitly shows convergence of the con- 
straints to zero, while Figure displays the associated 
convergence factors for those simulations. The latter are 
obtained by dividing the energy for different resolutions 
in pairs (i.e., A/'(Ai)/A/'(A2) and A/'(A2)/A/'(A3)). For a 
second order accurate code this convergence factor should 



be two in the convergent regime. FigureOshows that sec- 
ond order convergence is lost after some time, but this is 
expected for such long simulations, owing to accumula- 
tion of truncation error. However, the convergence factor 
gets closer to two when computed with the two highest 
resolutions. 



-10, 



— A 


= 0.025 


-- A 


= 0.0125 


... A 


= 0.00625 




,V,.'^1m'',I\i1|I III II ]- 



250 



time/(crossing-time) 



1000 



1250 



FIG. 1: This figure shows the logarithm of the energy for the 
constraints, in simulations of the gauge wave with amplitude 
A = 0.5. Three different resolutions are used, the coarsest 
run exhibits a slow growth which is negligible at higher reso- 
lutions. 
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FIG. 2: The logarithm of the relative error in the maximum 
value attained by gxx (compared to its exact value) versus 
number of crossing times for the simulations of Figure Q 



2. Amplitude A = 1 

Increasing the amplitude of the gauge wave introduces 
some complications when comparing to lower amplitude 
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FIG. 3: Convergence of the constraints energy to zero as a 
function of the number of crossing times, for the simulations 
of Fig. Inland Fig.Q This figure shows the convergence factors 
obtained by dividing the energy found at different resolutions 
in consecutive pairs. The oscillations observed are a product 
of phase velocities differing at different resolutions. 



runs at the same resolutions. For example, at the coars- 
est resolution (A — 0.025) the fields and the energy grow 
considerably: after 750 crossing times the energy is of 
order 10^. Clearly, even in the convergent regime, errors 
of this magnitude mean the numerical solution is of lit- 
tle use. Simulations with errors below ten percent last 
until about 600 crossing times for this particular grid 
resolution. To demonstrate the effect of resolution on 
the quality of the solution, we take 750 crossing times 
as the end-point of our simulations. Figures 01 El and 
El illustrate the results. The energy of the constraints, 
as shown in Figure^ shows marked exponential growth, 
with a large growth rate for the coarsest resolution. How- 
ever, increasing the resolution diminishes the growth rate 
considerably, and simulations can be extended for at least 
1500 crossing times (at which point we simply stopped 
the simulations, with small errors). 

Figure El shows the relative error in the maximum at- 
tained by gxx for the high amplitude gauge wave A = 1. 
Again, a rapid raise in the error is clearly seen for the 
coarsest resolution, but this effect is less noticeable at 
higher resolutions. At 750 crossing times, the error in 
the maximum of the function is less than ten percent. 
Finally, the convergence of the code is explicitly illus- 
trated in Fig. El which shows the convergence factors 
obtained by taking the energy for the constraints at dif- 
ferent resolutions and dividing them in pairs. The order 
of convergence of the constraints to zero is close to two 
for some time, and the lenght of this time increases with 
resolution. Finally, to illustrate the overall quality of 
the obtained numerical solution for different resolutions, 
figures [71 El and El display explicit comparisons between 
the exact solution and the middle and fine resolutions 
{p = 2,4). Figure [3 presents snapshots of the solution at 



280 and 560 crossing times, clearly both drifts in phase 
and amplitudes are observed for the middle resolution 
while for the finer one the difference is mainly observed 
in the phase. Despite these differences, the solution ob- 
tained is converging to the exact one. Figure El provides 
the convergence factor calculated with the two resolu- 
tions. As before, the factor deviates from its expected 
value after a while, but this improves with resolution. 
Finally figure E| shows the L2 norm of the difference be- 
tween the exact and numerical solutions. Clearly even 
after 600 crossing times, the fine resolution stays close 
to the exact solution and does not exhibit growth faster 
than the expected linear one. 
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FIG. 4: Logarithm of the energy for the constraints for the 
A — 1 gauge wave simulations. The coarsest resolution ex- 
hibits a marked growth after 600 crossing times, though with 
higher resolutions this effect becomes negligible. 
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FIG. 5: The logarithm of the relative error in the maximum 
attained by g^x for the simulations of Fig. |11 The errors 
exhibit a slow growth which diminishes with resolution. 
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FIG. 6: Convergence of the constraint energy to zero, for 
the simulations of Figures |1| and |3 obtained by dividing the 
energy for each resolution in consecutive pairs. As in Figure 
13 after some time the convergence factors deviate from the 
expected value of two, but this difference diminishes when the 
two highest resolutions are used to compute the convergence 
factor. 
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FIG. 7: Snapshots at 280 and 560 crossing times displaying 
the exact solution and the numerically calculated ones for the 
middle (A = 0.0125) and fine (A = 0.00625) resolutions. 
Convergence to the exact solution is evident as resolution is 
increased. 



A. Observations 

As mentioned above, the constraints C, Q, and C'ijki 
are initially satisfied to the level of round-ofF in the gauge 
wave tests. For grids with at least = 161 points, these 
constraints remain quite small throughout the runs, even 
for the high amplitude case, A = 1. This has several 
consequences: 

• First, most hyperbolic formulations differ, among 
other things, in how the constraints are added to 



FIG. 8: Convergent factor calculated with the numerical so- 
lutions for the middle and fine resolutions. At late times 
acumulation of errors makes the obtained value deviate from 
its expected one. 
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FIG. 9: L2 norm of error in the numerical solution. Although 
a faster than linear growth is observed for the middle resolu- 
tion, this effect converges away as a finer resolution is consid- 
ered. 



the right hand side of the equations 14] . Since the 
constraints themselves stay negligible small, we ex- 
pect the conclusions found with this test should 
be applicable to most hyperbolic formulations with 
the same choice of lapse and shift. Namely, that the 
use of symmetric hyperbolic systems and numerical 
techniques guaranteeing stability at the linear level, 
plus the addition of a small amount of dissipation 
stabilizes the problem. For instance, without the 
use of dissipation, for the case A — 1 and — 161 
points we can follow the system for about ten cross- 
ing times before the errors in the numerical solution 
become of order one. With a small dissipative term, 
on the other hand, the system after 1200 crossing 
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times does not yet exhibit errors of order one. 

• Second, the minimization technique presented 
in Section IIIII expUcitly makes use of non- 
homogeneous terms (in the free constraint- 
functions) in the expression for the energy growth 
in order to minimize it. However, in the current 
case these non-homogeneous terms are consider- 
ably smaller than the homogeneous contribution. 
Therefore, one would need to use huge constraint- 
function values for them to play a role in mini- 
mizing the growth, which would require the use of 
an extremely small Courant factor. Note that the 
formulation used employs the addition of the con- 
straints C, Ci and Ciju but not djk and = 
diN/N — Ai. The latter ones are not satisfied to 
round-off level but just to truncation level initially. 
As a result one could have considered a constraint 
energy which includes them or a modification of 
the formulation employed by the addition of these 
constraints in a suitable manner. We have not ex- 
plored these options in the present work as in the 
black hole case the non-homogeneous terms in the 
constraint energy are not negligible as in the gauge 
wave case. Indeed, this is what would be expected 
for generic scenarios. 



VI. BLACK HOLE SIMULATIONS 

We now turn our attention to the Schwarzschild black 
hole spacetime. 

In this case discretization errors (that is, due to fi- 
nite grid spacing) make the constraints start-off at non- 
negligible values and so the minimization technique can 
be used to pick up preferred formulations for this prob- 
lem. Before applying the minimization procedure we ex- 
amine the system with fixed constraint-functions with 
the goal of understanding some specific issues. As in 
the rest of the paper, the simulations of this section are 
done with the two-constraint-function family of formula- 
tions discussed in Section II (the reason for not using the 
single-constraint-function family is explained in Section 
IVI A 1|) : and the fixed values for 7 and rj here used are, 
in the absence of any other obvious choices, 7 = = 77. 
Recall that this two-constraint-function family is sym- 
metric hyperbolic for any values of rj and 7, in particular 
for 7 = = ?7 

Next we study the infiuence of the position of the in- 
ner boundary; the outer boundaries (their influences be- 
ing discussed later in the paper) are kept at the same 
position, ±5M in all the runs of this section. The outer 
boundaries are chosen quite close on purpose, to make 
the turn-around time for the runs shorter for the rather 
detailed analysis of this section. However, it should be 
clear that the points here made are quite independent on 
the particular position for the outer boundaries chosen. 
Later, when applying the minimization of the constraints, 



we will choose several different values for the position of 
the outer boundary. 



1. Inner boundary and the outflow condition 

Black hole excision is usually based on the assump- 
tion that an inner boundary (IB) can be placed on the 
domain such that information from this boundary does 
not enter the computational domain. The boundary is 
supposed to be contained inside the black hole and also 
to be purely outflow, i.e., all modes propagate off of the 
grid at the boundary. This requirement places strenuous 
demands on cubical excision for a Schwarzschild black 
hole in Kerr-Schild, Painlevee-GuUstrand or the Martel- 
Poisson j22| coordinates: the cube must be inside 0.37M 
in each direction. This forces one to excise very close to 
the singularity, where gradients in the solution can be- 
come very large, requiring very high resolution near the 
excision boundary to adequately resolve the solution. Fi- 
nally, we note that this requirement follows directly from 
the physical properties of the Schwarzschild solution in 
these coordinates, and is independent of the particular 
formulation of the Einstein equations 0. 

With our current uniform Cartesian code, we are not 
able to provide the resolution required to adequately rep- 
resent the Schwarzschild solution close to the singularity. 
While we are actively working on solutions to this prob- 
lem, currently our only practical alternative is to place 
the inner boundary inside the event horizon, but outside 
the region specified by the outflow condition, i.e., the so- 
lution has incoming modes on the inner boundary. One 
could attempt to provide data for the incoming modes on 
the excision boundary. However, there is no general the- 
ory well-posed problems when the principal part rank, 
i.e., number, of zero speed modes on the boundary, is 
not constant. Moreover, ill-posed problems for such con- 
figurations are known. Thus we simply do not apply 
boundary conditions to the incoming modes on the inner 
boundary, resulting formally in an ill-posed problem . 
In this section, however, we argue that errors from this in- 
consistency do not prevent us from learning much about 
the numerical properties of our formulation for black hole 
spacetimes. 

Figure H10|l shows the results of simulations with differ- 
ent positions for the inner boundary (IB), with A = M/5, 
obtained without imposing boundary conditions at the 
IB. The IB at 0.3A/ (i.e. hah the length of a cube 
centered at the origin) corresponds to a purely outflow 
boundary. At the other extreme, the inner boundary at 
1.3M gives an inner boundary that penetrates outside 
the event horizon. The cases between 0.3A/ and 1.3M 
correspond to IB inside the black hole, with inflow por- 
tions. All runs shown in Fig. ^1 except the first with 
IB at 0.3M, should have convergence problems, since not 
giving boundary conditions is inconsistent with the struc- 
ture of the characteristic modes. However, one also sees 
that at this resolution, placing the inner boundary so 
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close to the singularity causes the code to run a factor of 
ten less, presumably because of lack of resolution. 

Different positions for the inner boundary 



Convergence test, fixed constraint-functions {y=0. r\=0) 
Energy for the constraints 




100 



FIG. 10: Black hole simulations, with fixed outer boundaries 
(at 5M), and different position for the inner boundary (IB). 
The resolution is A = A//5. 



Figure (|ll|l shows a convergence test for a configura- 
tion test with a the excision boundary at I.IM with reso- 
lutions A = M/5, M/10, M/20. From here on the errors 
shown in the plots are those of the numerical solution 
Un relative and with respect to the exact one ite; more 
precisely, the L2 norm of 



-1/2 



where the sum is over the components of the vector val- 
ued functions u. 

This test can give some indication of possible numer- 
ical stability problems. While the solution diverges in 
all cases, the fact that the code appears to converge in 
the short term indicates that at these resolutions, and 
for these run times, the expected instability owing to im- 
proper boundary placement appears to grow slower than 
other unstable modes in the solution. Thus we can still 
obtain valuable information about the solution and its 
numerical properties. Finally, we emphasize that the 
inconsistent inner boundary probably leads to conver- 
gence difficulties that could be detected with more exten- 
sive tests, such as performing a Fourier decomposition of 
the numerical solution, and a convergence test frequency 
by frequency likely would make the numerical instability 
manifest, see for instance Q- 



A. Dynamic minimization: preliminary discussion 

In this section we examine several issues that arise in 
the constraint energy minimization procedure. 
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FIG. 11: Two-constraint-function family, with fixed values 
7 = = 77, inner and outer boundaries at l.lAf and 5AI, 
respectively. 



1. Why one should use the two- constraint- function family 

We consider only the two-constraint-function family 
of formulations of the equations for all constraint min- 
imization runs. The single-constraint-function formula- 
tion is inadequate because, as mentioned above, the free 
constraint-function x may not equal zero. Thus x may 
only be negative or only positive during an entire run to 
be a continuous function of time. This considerably limits 
the power of the dynamic minimization technique, since 
in order to control the constraints one might need at some 
time a positive value of x some other time a neg- 

ative one. Indeed, this occurs as Figure H12|l shows. For 
this resolution and location of inner and outer boundaries 
the initial discretization error for the constraints leads to 
a value of the energy of 



7V(0) 0.99925 x 10" 



(33) 



Figure [T^ thus shows two evolutions, one with a negative 
seed value, x — —1-0, and a tolerance value 10~^. The 
second run has a tolerance value of 10~*, and as seed 
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value X = 1.0. In both cases Ha = 1- In both cases 
X changes sign, indicating that this can be expected in 
general, and that there is no continuous interpolation for 
x{t) in the limit A = such that x 7^ 0- 

Dinamio % and its sign change 




FIG. 12: Generally x would need to change sign in order to 
achieve certain control on the constraints. Therefore, it would 
not have a continuous limit when At 0. For this reason, 
the two-constraint-function formulation is used throughout 
this paper. 



the spatial derivatives are discrete, but time continuous. 
While a fully discretized method could be developed, we 
simply use the semi-discrete analysis here. In the limit 
Ai — > 0, one naturally expects this semi-discrete analysis 
to be a perfect description of the fully discrete scheme. 
Here we verify that the fully discrete evolutions are, in- 
deed, very well approximated by the semi-discrete anal- 
ysis, cf. Eq. H12|) . even for rather large Courant factors, 
such as A — 0.5 that we use for the black hole runs of 
this paper. 

Figure 1131 shows results of an evolution with inner 
boundary at I.IA/, A = 0.5,(7 = 0.03, outer boundaries 
at ±5M, and A — M/5. There are two curves for the 
time derivative of the energy, JV. The first curve is ob- 
tained via the semi-discrete prediction, Eq. (|12l) . with the 
matrix valued integrals computed during evo- 

lution. The second curve is obtained from second order, 
centered differences of the energy, Af. Both curves agree 
remarkably well, which indicates that the semi-discrete 
expression for Af captures the dynamic constraint behav- 
ior in the fully-discrete simulation with a Courant factor 
of A = 0.5. 

As another illustration of this, figure E| shows the 
results obtained for one of the gauge wave tests presented 
above: A = 0.5, 161 points, and A = 0.25. For a clean 
visualization of the agreement we show here a period of 
time between 30 and 50 crossing times. The agreement 
of the curves is evident. 



2. The accuracy of a semi-discrete picture. 
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FIG. 13: Energy for the constraints and its time deriva- 
tive, J\f, computed through the semi-discrete prediction and 
through numerical differentiation. The remarkable agreement 
between the two indicates that the semi-discrete analysis used 
to calculate the constraint minimization is faithful represen- 
tation of the fully discrete evolution. 



The constraint minimization method, as described in 
Section IIIII is based on semi- discrete equations, where 
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FIG. 14: Time derivative of the energy of the constraints 
computed through the semi-discrete prediction and through 
numerical differentiation for the gauge wave case. The agree- 
ment of the two curves is evident. 



3. Practical questions regarding constraint minimization 

We now consider some practical questions that arise 
when performing the constraint minimization method, 
namely: how often should the minimization be per- 
formed?, and how fast should the constraint-functions 
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FIG. 15: The plots in this figure and those in Fig. 1161 show 
the eflects of varying the frequency of performing the con- 
straint minimization, as well as the dependence on Ua- The 
constraint-function rj for na — 10* is not shown as the run 
crashes very early and the scale in the time axis for the asso- 
ciated plot is not logarithmic. 



multiplying the constraints in the RHS be allowed to 
change? The minimization procedure can be computa- 
tionally expensive, making it advantageous to perform 
the analysis infrequently if possible. (Note that the eval- 
uation of the constraint energy require knowledge of the 
right hand side of the equations, so each evaluation takes 
about as much computer time as is required in a time- 
step). 

The second question relates to how fast the constraint- 
functions are allowed to vary by setting , which is used 
to estimate how many steps are required for J\f to relax to 
the tolerance value, T. Furthermore, can one trade the 
frequency of the constraint minimization with different 
values of na? For example, is performing the minimiza- 
tion at every iteration with a large value of Ua equivalent 
to performing the minimization every certain number of 
steps with a smaller value of Ual. We will see that it 
is preferable to perform the minimization at every time 
step. 



Figure (|15|l shows the results of runs with the same 
numerical constraint-functions as in the previous subsec- 
tion: outer boundaries at ±5M, inner boundary at I.IM, 
a = 0.03, A = 0.5, 51'^ points, and a tolerance value 
T ~ 10^^, performing the minimization at every time 
step, but now with different values of na- The constraint- 
function 7 is fixed to 7 = 0, and the minimization of the 
constraints is applied using rj, as described in Section lTTll 
The constraint energy, Af is shown in the upper panel, 
and r]{t) in the lower panel. With a fixed Courant fac- 
tor, small values of Ua are problematic because large and 
fast variations in r]{t) are allowed and, as shown in the 
figure, do occur. On the other hand, large values of Ua 
can let the energy grow too much. Compare now Fig- 
ure (|16(l . where the minimization is now applied every 
ten time steps. As can be seen, it is better to apply the 
minimization at every time step. Otherwise the energy 
appears to grow too fast between recalculations of the 
constraint-functions. Finally, we note that this behavior 
may be model dependent, and in other scenarios it may 
be possible to use the minimization less frequently. In 
this paper, however, we perform the minimization after 
every time iteration 1361 . 



4- Sensitivity to the tolerance value 

We choose constraint-functions in the constraint min- 
imization procedure such that the energy, Af, decays to 
a tolerance value, T, after certain number of time steps, 
Ua- In this section we discuss reasonable choices for T, 
and discuss its the influence on the final solution. There 
are some reasons to believe that a value close to the ini- 
tial discretization error is a good choice 0|. Here we 
present additional evidence for this by comparing simu- 
lations with different values of T and Ua- 

Figures El and El show A/" as a function of time, for 
T = 10-2, 10-3, 10"^, 10"^ respectively, each one with 
Ua = 1, 10, 10^, 10^, 10"*. It can be seen that one can in- 
deed perform similarly choosing a tolerance value that 
is below the discretization error by using an appropriate 
value for n^. But notice that the longest runs are ob- 
tained when M naturally has a value near the initial dis- 
cretization error. Even if T is very small, a large value for 
na allows the constraints-functions to change only very 
slowly, resulting in slow changes of J\f towards T. In 
summary, a combination of T and Ua that keeps J\f near 
the initial discretization error appears to give the best 
results. 



5. Limitations of the constraint-minimization procedure 

Fig. (|15|l shows that the constraint minimization gives 
an improvement of five to ten times in the lifetime of the 
simulations, compared to the results shown in Fig. H10|l . 
the method does not prevent the eventual code crash 
when J\f is very large. Given that the constraint mini- 
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FIG. 16: This plot shows simulations as those of Fig. (I15L 
except that here the minimization is done every 10 iterations. 
Comparing it with the previous plot, it seems clear that per- 
forming the minimization at every iteration seems a better 
option 



FIG. 17: This Figure and Figure lTSl demonstrate the influence 
of the tolerance value, T, on the results of the constraint min- 
imization. The two-constraint-function formulation is used, 
with 7 = and r){t) chosen dynamically. The results are 
not sensitive to the chosen tolerance value, provided that one 
avoids large variations in ri{t) by using appropriate values of 
Ua- However, notice that "appropriate" values of Ua naturally 
keep the value of the energy close to its initial, discretiza- 
tion value (here given by 0.99925 x 10~^. Therefore, keeping 
the constraint energy near its initial discretization value ap- 
pears to give the best results. This figure shows runs with 
T = 10"^ and T = 10"^. Fig. IT^ shows runs with T = 10"* 
andT= 10"^ 



mization method is designed to prevent this, we consider 
possible reasons why the method eventually fails. Oper- 
ationally, the failure seems to occur because large scale 
variations of 77 on ever decreasing time scales are required 
at late times. It is not possible for the code to resolve 
such variations in the fields with a fixed Courant factor. 
This was partly analyzed in Section IVI A 31 where we 
discussed the dependence of ?]{t) on Ua (see Figure [T^. 
Furthermore, when is large, M does not return im- 
mediately to a value near T, but changes are affected 
slowly. Thus, if jV grows on time scales of a large num- 
ber of steps, the minimization method is unable to halt 
the growth. Having evidence of why the code crashes, 
we now can attempt improvements o the method, as dis- 
cussed in the next section. 



B. Two dimensional minimization and numerical 
results 

In this section we exploit the freedom in the two- 
constraint- function family to extend the lifetime of black 
hole simulations. The tolerance value for the energy is 
chosen to be a value roughly one order of magnitude 
larger than the initial discretization error, and Ha is set 
to either 10^ or lO'^. The boundaries are placed at 5M, 
as in the runs discussed previously, and also at lOM and 
15M. 

We exploit the fact that there are two free constraint- 
functions to achieve not only a given tolerance value, but 



15 



10r 

10° r 

^ -2' 
i 10 r 

LU 

10-^ r 
10" 







10 r 



10 r 



g 10 

LU 

10"' 



10' 



Tolerance=10" 
Fixed y, dynamic t\ 



Minimization every one thousand iterations 



n^=1 

n^=10= 
n =10" 



10 



time[M] 

Tolerance=10"^ 
Fixed y, dynamic t\ 



100 




1 10 
time[M] 



100 



1000 




40 
time [M] 

Minimization every one iteration 




FIG. 18: This figure compares results for T = 10 ^ and T = 
10~^. See Fig. llTl for additional information. 



also to minimize the change in the constraint-functions 
r]{t),"/(t), to prevent fast variations in these constraint- 
functions, as explained below. Thus, with two constraint- 
functions to achieve the desired tolerance value, we can 
impose an additional condition to minimize the varia- 
tion in the constraint-functions from one time step to 
the other one. The motivation for this condition comes 
from the discussion in the previous section on the limi- 
tations of the constraint minimization technique, where 
large oscillations in the constraint-functions were needed 
in order to keep the constraints under control. Therefore, 
it seems reasonable at this stage to conjecture that the 
lifetime would be extended even more if one was able to 
apply the constraint minimization in a way such that fast 
variations in the constraint-functions are not needed. As 
shown below, this conjecture seems to be correct. There- 
fore, within all the constraint-functions that achieve the 
desired energy growth for the constraints, we choose at 
time step n + I the pair that minimizes the quantity 

A := [-nin + 1) - Tj{nf + [7(71 + 1) - j{n)f (34) 

To apply this condition, consider that Eq. (|21|l indi- 
cates that Af is non-linear in 7 but linear in rj, allowing 



FIG. 19: These plots show the associated rj{t) for the Ua = 
1, 1000 cases of the previous plot, in more detail. Notice the 
scale in the ria = 1 case, rj changes a lot per time step. On 
the other hand, notice how in the Ua = case ri{t) changes 
slower but still changes, and eventually it needs to take very 
large values in order to control the constraints. 



one to solve for rj such that Af = —aAf, 

_ -{aAf + l''""' + 27) + 27TX 

X"(l + 27)+7Xx ^-^^^ 

where, as in Section III, a is given by Eq. ljlTII . 

A set of values for 7 is chosen within some arbitrary, 
large interval. For each 7 the corresponding t] given by 
Eq. (|35|l is computed, and the pair (77, 7) that minimizes 
A defined in Eq. H34|) is chosen. As explained in Section, 
the constraint-functions (7, 77) can take any value, except 
for 7 = —1/2, value for which the equations are singu- 
lar. Therefore, there are two "branches" , we have only 
explored the one associated with 7 < —1/2, by using as 
seed values 77 = 0,7 = —1 and restricting the minimiza- 
tion procedure to values 7 < — 1/2. As discussed next, 
almost an extra order of magnitude in the lifetime of the 
simulations can be obtained, and the run that initially 
lasted for lOM without the minimization of the con- 
straints, now runs for around 700M — lOOOAf (without 
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any symmetry 
imposed). 



bitant, octant, or of any other type- 
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FIG. 20: Dynamic minimization done with boundaries at 5M, 



A = M/5, T= 10"^ and Ua 



Figure 1201 shows the results of a simulation with res- 
olution A — M/5, T — 10~^ (close to the initial dis- 
cretization error value, given by Eq. OSI) and Ua = 10^. 
The figure shows that as the code crashes the constraint- 
functions start having large and fast variations. A nat- 
ural question that this raises is whether these varia- 
tions are a cause or consequence of the code crashing. 
For reasons discussed below, they appear to be a conse- 
quence. Figure UT\ shows the same run as that shown in 
Fig. 1201 except that the minimization is stopped at 750M 
(at which point the constraint-functions are, nearly con- 
stant), and from there on the last value of the constraint- 
functions is used; namely. 



f] = -1.88,7 : 



-1.00 



(36) 



The code still crashes at roughly the same time. There- 
fore, the variations in the constraint-functions observed 



in Fig. 1201 do not cause the code to crash, but appear 
to be a consequence of other instabilities. Figure |221 
shows a convergence test with two resolutions (A — M/5, 
A = M/10), and keeping the constraint-functions (ob- 
tained from the A = M/5 resolution run) fixed after 
750M. 

Dynamic minimization, boundaries at 5M 
Keeping tfie constraint-functions fixed after U750M 



C. functions l<ept fixed after 750M 
C. functions not kept fixed 
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FIG. 21: Same as previous Figure, but keeping the constraint- 
functions constant after 750M. The figure compares the re- 
sulting energy for the constraints with that of the previous 
figure (shown at late times only, since because of the setup 
the runs are identical up to t = 750M). 



Another measure of the error in the solution is the mass 
of the apparent horizon. Figure 1231 shows this mass us- 
ing the dynamic constraint-functions obtained from the 
simulation of Figure 1201 and by running at each iteration 
Thornburg's apparent horizon finder j24| . For the coars- 
est resolution, the initial value of the mass, as given by 
the horizon finder, is l.OOTAf. Compared to this value, 
the initial oscillations have a relative error of less than 
one percent. After some time, the mass approximately 
settles down to a value that is around 1.009, which cor- 
responds to an error of the order of one part in one thou- 
sand. For the higher resolution, the initial value of the 
mass as given by the horizon finder is 0.99951. With re- 
spect to this value the initial oscillations are at most of 
the order of one part in one thousand, and at late times 
the apparent horizon mass settles down to 0.99953, which 
corresponds to a relative error of one part in 10^. 

Even though the constraint-functions do not settle 
down completely to a stationary value, they oscillate very 
little. One question that this raises is how does the code 
perform if one fixes these values given by Eq. (|36|l from 
the very beginning. The plots in figure 1241 make this 
comparison. Interestingly, the run with dynamSic mini- 
mization runs slightly longer, even when the constraint- 
functions after some time are essentially constant, and 
even though the solution being evolved is stationary at 
the analytical level. This shows that the dynamic mini- 
mization not only requires little experimentation but also 
seems to be effective in that it naturally allows for varia- 
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FIG. 22; Convergence test, with dynamic constraint-functions 
kept constant after 750A/. 



tions in time that accommodate to the variations of the 
numerical solution. 



2. Boundaries at lOM. 
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FIG. 23: Apparent horizon mass, with dynamic constraint- 
function values. The run with higher resolution ran out of 
computing time (in the equivalent plots of the previous Figure 
an apparent horizon was not searched for). 



from the run with dynamic minimization and boundaries 
at 5A/ do not perform as well as those obtained with 
boundaries at lOAf . However, even using the constraint- 
functions obtained from the simulation with boundaries 
at 5Af is much better than using a naive choice (say, 
7 = = 0, which for the resolution of Fig. 1261 runs for 
less than 30M, as shown in Fig.lllH. 

Figure 1271 shows the apparent horizon mass, for the 
simulations of Figures 1251 and one simulation from 
Fig. 1261 with constraint-function values given by Ea. (|37|l . 
In both cases the resolution is coarse, A = A//5. As 
for the case with boundaries at 5Af and with the same 
resolution, the errors are less than one percent when com- 
pared to the mass given by the initial data. From Fig. [57| 
one can also see that the oscillations in the mass do not 
seem to be caused by the time variation of the constraint- 
functions, as they are still present in the case in which 
fixed constraint-functions arc used. 



We now examine data from a configuration equivalent 
to that of the previous section, except that now that the 
boundaries are at lOA/. The initial discretization value 
for the energy is N[Q) = 1.2845 x 10-^ and T = lO"", 
and Ua — 10^ are chosen. As seen in the previous case, 
the constraint- functions eventually settle into oscillations 
about fixed values, 

ri = -2.96 X 10"\ 7 = -2.48, (37) 

as shown in Figure 123 These steady-state values are 
quite different from the previous configuration with 
boundaries at 5M, see Ea. (|5^ . This raises the question 
of what would happen if one ran with boundaries at 10 A/, 
and fixed constraint-functions given by Ea. (|37() in one 
case and Ea. (|36|l in the other. Figure I^Hl makes this com- 
parison. As expected, the constraint-functions obtained 



3. Boundaries at 15M . 

Finally, we consider a configuration with boundaries 
at 15M, though with less detail as before. Figure [SS] 
shows results data equivalent to those discussed for Fig- 
ures |2D1 and I^Hl except now that the boundaries are at 
15Af. The initial, discretization value for the energy is 
7.6459 X 10-6, and T = 10-^ = 100 was used. The 
minimization of the constraint-functions is stopped at 
450 Af, at which point the constraint-functions are ap- 
proximately constant, and equal to 

,y = -1.35 X 10^1 , 7 = -3.39. (38) 

Figure |2H1 shows that the dependence of the lifetime 
on the location of the outer boundaries is not monotonic, 
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FIG. 24: Comparison between dynamic and fixed, fine-tuned 
constraint-function values given by the asymptotic values at 
which they approach [cf. Eq. (1361 1. 



FIG. 25: Dynamic minimization done with boundaries at 
lOM, A = M/5, T = 10"", and Ua = 10^ 



as for this case the code runs for, roughly, lOOOM, while 
with boundaries at lOAf and 5M it ran for around 700M, 
and 800A/, respectively. A detailed analysis of such de- 
pendence would be computationally expensive and be- 
yond the scope of this work, and may even depend on 
the details of the constraint minimization, such as the 
values for T and Ua- 



VII. FINAL COMMENTS 

This paper presents a number of new techniques ap- 
plied to the simulation of Einstein equations, namely: (1) 
a symmetric hyperbolic formulation with live gauges; (2) 
a numerical discretization based on the energy method 
with difference operators that satisfy summation by 
parts and a projection method to apply boundary condi- 
tions; (3) a constraint minimization method for dynami- 
cally choosing constraint-functions that multiply the con- 
straints in the evolution equations without requiring spe- 
cial knowledge of the solution. j37| 

We use a generalization of the Bona-Masso slicing con- 



ditions, and to date, this is the only formulation of Ein- 
stein's equations with this slicing conditions known to 
be symmetric hyperbolic (for symmetric hyperbolic for- 
mulations with other dynamical gauge conditions see 
|2It|V There are strongly h yper bolic formulations with 
this gauge (see, for example, l26| and references therein), 
though, recall that contrary to some common belief, 
strong hyperbolicity does not automatically define a well- 
posed initial value problem (IVP). A well-posed IVP for 
strongly hyperbolic can be found by requiring the exis- 
tence of a smooth symmetrizer. However, this smooth- 
ness is a non-trivial condition and it is usually not stud- 
ied in formulations of Einstein's equations. Some al- 
gebraic conditions do imply the existence of a smooth- 
symmetrizer |23|, but for the Bona-Masso slicings these 
conditions can only be a priori guaranteed for the time- 
harmonic subcase |l3| . In the presence of boundaries the 
situation is even more complicated: there are examples in 
the context of Einstein's equations explicitly showing ill 
posedness of certain strongly hyperbolic equations which 
do have smooth symmetrizers, when maximally dissipa- 
tive boundary conditions are used (while for symmetric 
systems such a problem is known to be well posed) jg^. 
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FIG. 26: Running with boundaries at lOM, using fixed, 
but fine tuned constraint-functions, obtained from runs witfi 
boundaries at 5M and lOM. 
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FIG. 28: Dynamic minimization done witli boundaries at 
15M, A = M/5, T = 10"^ and Ua = 10^ Tlie constraint- 
functions are constant for t > 450M, where they are rj = 
—0.135,7 = —3.389; thus, there is no response of the con- 
straint functions when the code is about to crash. 



FIG. 27: Apparent horizon mass for the simulation of Figure 
1251 and the simulation of Fig. 1261 with constraint-functions 
given by Eg. 1371 . A logarithmic scale in time is used in or- 
der to show the oscillations in more detail. The oscillations 
are not caused by time variations in the constraint-functions, 
since they are also present in the fixed-constraint-functions 
case. 



While we use time harmonic slicing in this paper, the 
freedom to use other slicing conditions could prove use- 
ful in other scenarios [29l |. 

Finite difference discretizations based on the energy 
method A, J5, exploit results that rigorously guarantee 
linear numerical stability of IVPs as well as IB VPs. In 
particular, stable simulations of the gauge wave with pe- 
riodic boundary conditions are obtained for large values 
of amplitude, A, for at least a thousand crossing times. 
These simulations show that the constraints remain well 
behaved throughout the evolution, indicating that con- 
straint violations, if any, grow very slowly in time. For 
the black hole cases, where both inner and outer bound- 



aries are present, these numerical techniques allow for 
a clean handling of particularly delicate issues. For in- 
stance, defining the difference and dissipative operators 
at these boundaries, and how boundary conditions are 
imposed (in particular, in non-smooth boundaries) are 
addressed in a completely systematic way. 

Numerical stability guarantees that errors do go away 
with resolution, but at fixed resolution they can still grow 
quickly in time. These fast growing errors can be intro- 
duced by the continuum instabilities, triggered by numer- 
ical errors, by the numerical scheme, or an y co mbination 
of the above. The technique explored here [ij, automat- 
ically adjusts the formulation of Einstein's equations in 
such a way that the discrete constraint violations follow 
some prescribed behavior (for example, their norm re- 
mains close to its initial, discretization value). There are 
a number of lessons learned from the application of this 
dynamic minimization procedure in this paper, which are 
worth highlighting: 

1. The semi-discrete picture describes what happens 
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in the fully discrete case remarkably well, even for 
cases that are highly non-stationary. Here by semi- 
discrete one means a picture that assumes time to 
be continuous, but space to be discretized with an 
arbitrary (not necessarily small) grid spacing. 

2. The technique of can be used not only as a 
practical tool for extending the lifetime of the sim- 
ulations, but also for gaining conceptual insight in 
the problem of constraints violations in free evo- 
lutions. Sometimes very large adjustments must 
occur on short time scales in order to control the 
constraints, which may become too fast or large for 
a fixed Courant factor. There are two issues related 
to this observation: 

(a) Very likely this same feature is present in 
many other cases, where different formulations 
of the equations and numerical techniques are 
used. It clearly points out a limitation in ad- 
justing the equations so as to minimize the 
constraints growth, independently of what the 
adjustment technique is. 

(b) Nevertheless, the identification of this limita- 
tion points out a way to proceed with the tech- 
nique of j^. Namely, to take advantage of 
many-constraint-function formulations by re- 
defining the equations in a way such that not 
only a given behavior for the constraints is 
achieved, but also the adjustment varies as lit- 
tle as possible between two time steps, as done 
in Section rvTBl 

The results of Section IVl Bl confirms to a good ex- 
tent the validity of the previous discussion. 

As a practical matter the lifetime of the full 3D black 
hole simulations are extended from around 20M up to 
700M— lOOOAf. This is achieved without employing sym- 
metry restrictions (like octant, bitant or any other), or 
previous knowledge of the expected solution. When cm- 
ploying symmetries the code actually runs much longer. 
For example, Figure shows a fully 3D simulation with 
boundaries at lOM, using the values around which the 
parameters settle down after a while, given by Ea. H37|) 
(that is, one of the simulations of Fig l26|l . compared to 
the same simulation in octant symmetry, which for con- 
venience was stopped at 18,000Af. 

At this point it is not clear why the code crashes af- 
ter 700M - lOOOM. One possibihty could be the pres- 
ence of a numerical instability caused by the inconsis- 
tency of having an inner boundary that is not purely out- 
flow, but is treated as if it was. Although we have done 
some convergence tests, numerical instabilities are some- 
times subtle, and very detailed and careful convergence 
tests have to be done in order to detect them, especially 
when the initial data has few and low frequencies (as is 
here the case) However, an important feature of all 
the simulations here presented, including those of Section 
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FIG. 29: This Figure shows a fully 3-D run with fixed and 
fined-tuned parameters, given by Eg. 13711 . compared to the 
same run in octant symmetry, which was stopped at 18, OOOM. 



WIBl is that large values of Ua had to be used in order 
to prevent large and quick variations in the constraint- 
functions. Therefore, one is not completely controlling 
the constraints — for that a value of tIq of order one would 
have to be used — and they do grow. Therefore, one pos- 
sibility for achieving small values of Ua without large and 
quick variations in the constraint-functions would be to 
introduce more free constraint-functions and to make use 
of this extra freedom as in Section lyiBl The results of 
Section lVl BI stronglv suggest that this should extend the 
lifetime even more, but more work must be done in or- 
der to explicitly study this. Finally, the constraint mini- 
mization method is designed to dynamically control the 
constraints growth without any a priori knowledge of the 
solution. Therefore, a natural next step also seems to be 
an application to dynamical spacetimes. 
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[35] There is a slight abuse of notation here, in the sense 
that a does not denote an index, as before. Similarly, the 
subscript in Ua indicates that the quantity is related to 
a through Ea. l|17|l . 

[36] At every full time iteration. That is, we do not perform 
the minimization at the intermediate timesteps of the 
Runge-Kutta integration, we have not explored this pos- 
sibility. 

[37] Some alternatives to control the constraints include ex- 
plicitly solving them instead of some of the evolution 



equations {constrained evolution, see for example [30(1 for 
recent work), projecting the solution onto the constraint 
surface {constraint projection [sj), enlarging the system 
in a way so as to force the solution to decay towards 
the constraint surface {lambda systems [s^), or adopting 
appropriate discretizations and/or gauge choices so that 
the discrete constraints are exactly preserved during evo- 
lution [13. 



